

 
% warning off
optz = optimset('fminunc');
optz = optimset(optz,'Display','off'); 
optz = optimset(optz,'Largescale','off'); 



%[aa,bb,cc] = tblreadrjv('cloilfutleducf2019.csv','comma');
%save cloilfutleducf2019 aa bb cc 
 load cloilfutleducf2019 aa bb cc 
 
 
 dcc = datenum(cc);
firstyear = min(year(dcc));
finalyear = max(year(dcc));
finalmonth = month(dcc(end));

fut1dtz = dcc;
bb(10,:)
fut1daily = log(aa(:,10));
ftfrntdaily = log(aa(:,1));
ftsprd = fut1daily-ftfrntdaily ;
[mftsprd ,mftsprdnan,mfutsprddate] = bettermonthavg([ftfrntdaily ftsprd],fut1dtz );

estdata = [diff(mftsprdnan(:,1)) mftsprdnan(2:end,2)];
estdates = mfutsprddate(2:end,:);
 ix = find(isnan(estdata(:,2)),1,'last');
 estdata = estdata(ix+1:end,:);
 estdates =estdates(ix+1:end);
 datestr(estdates(1,:))
 
 jmppoint = mod(3-month(estdates(1)),3); %to make sure that the first obs is end of quarter. 
 

 qestdata3 = estdata((jmppoint+1):3:end,:);
 qestdates =  estdates((jmppoint+1):3:end,:);

if mod( month( qestdates(1,:)),3) ~= 0;
    disp('Got the start date wrong')
end


stpoint = find( qestdates == datenum(2001,3,1));
enpoint = length(qestdates);
svvalues = NaN(enpoint,5);
 
for zx = stpoint:enpoint;
    newsv = [0.0032  0.7078    0.1049    0.1328 0.1];
    rvll =   @(params)llf_restat(qestdata3(1:zx,:)' ,params,[],[],4);
    outvalues = fminunc(rvll,newsv , optz);
    svvalues(zx,:) =    outvalues;
    
  
 
    
end

sv2values = NaN(enpoint,5);
for zx = stpoint:enpoint;
  rvll2 =   @(params)llf_restat(qestdata3(1:zx,:)' ,[0 params],[],[],4);
    outvalues = fminunc(rvll2,newsv(2:end) , optz);
    sv2values(zx,:) =   [0  outvalues];
end



sv3values = NaN(enpoint,5);
for zx = stpoint:enpoint;
  rvll3 =   @(params)llf_restat(qestdata3(1:zx,:)' ,[0 params 0],[],[],4);
    outvalues = fminunc(rvll3,newsv(2:end-1) , optz);
    sv3values(zx,:) =   [0  outvalues 0];
end


for zx = enpoint:-1:stpoint;
    
    rvorig =   @(params)llfpermtempdecompdiff(qestdata3(1:zx,1)',params);
    rvbout = fminunc(rvorig ,newsv(2:end-1),optz);
    sv4values(zx,:) = rvbout ;
    
end




figure
subplot(3,2,1)
plot(qestdates,svvalues(:,1))
title('Risk Premium Persistence')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,2)
plot(qestdates,svvalues(:,2))
hold on
plot(qestdates,sv2values(:,2),'g-')
hold on
plot(qestdates,sv3values(:,2),'r*-')

plot(qestdates(stpoint:enpoint),sv4values(stpoint:enpoint,1),'kp-')

title('Temp Shock Persistence')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,3)
plot(qestdates,abs(svvalues(:,3)))
title('Perm Shock Std. Dev')
hold on
plot(qestdates,abs(sv2values(:,3)),'g-')
hold on
plot(qestdates,abs(sv3values(:,3)),'r*-')
plot(qestdates(stpoint:enpoint),abs(sv4values(stpoint:enpoint,2)),'kp-')

datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,4)
plot(qestdates,abs(svvalues(:,4)))
title('Temp Shock Std. Dev')
hold on
plot(qestdates,abs(sv2values(:,4)),'g-')
hold on
plot(qestdates,abs(sv3values(:,4)),'r*-')
plot(qestdates(stpoint:enpoint),abs(sv4values(stpoint:enpoint,3)),'kp-')

datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,5)
plot(qestdates,abs(svvalues(:,5)))
hold on
plot(qestdates,abs(sv2values(:,5)),'g-')

title('Risk Prem Shock Std. Dev')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,6)
plot(qestdates,NaN+abs(svvalues(:,4)))
title('Legend')
hold on
plot(qestdates,NaN+abs(sv2values(:,4)),'g-')
hold on
plot(qestdates,NaN+abs(sv3values(:,4)),'r*-')
plot(qestdates(stpoint:enpoint),NaN+abs(sv4values(stpoint:enpoint,3)),'kp-')
legendhandle =legend('New Code Rho_risk unconstrained','New Code Rho risk zero','New Code Rho risk = sig risk = 0','Old Code Single Variable')
set(legendhandle,'fontsize',12)
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


orient tall
print -dpdf lotsofstuff




figure
subplot(3,2,1)
plot(qestdates,svvalues(:,1))
title('Risk Premium Persistence')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,2)
plot(qestdates,svvalues(:,2))
 title('Temp Shock Persistence')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,3)
plot(qestdates,abs(svvalues(:,3)))
title('Perm Shock Std. Dev')
 datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,4)
plot(qestdates,abs(svvalues(:,4)))
title('Temp Shock Std. Dev')
 datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,5)
plot(qestdates,abs(svvalues(:,5)))
 title('Risk Prem Shock Std. Dev')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])




orient tall
print -dpdf simplestuff



%Try Again in Dollars!


clear
load cloilfutleducf2019 aa bb cc 
 optz = optimset('fminunc');
optz = optimset(optz,'Display','off'); 
optz = optimset(optz,'Largescale','off'); 


 
 dcc = datenum(cc);
firstyear = min(year(dcc));
finalyear = max(year(dcc));
finalmonth = month(dcc(end));

fut1dtz = dcc;
bb(10,:)
fut1daily =  (aa(:,10));
ftfrntdaily =  (aa(:,1));
ftsprd = fut1daily-ftfrntdaily ;
[mftsprd ,mftsprdnan,mfutsprddate] = bettermonthavg([ftfrntdaily ftsprd],fut1dtz );

estdata = [diff(mftsprdnan(:,1)) mftsprdnan(2:end,2)];
estdates = mfutsprddate(2:end,:);
 ix = find(isnan(estdata(:,2)),1,'last');
 estdata = estdata(ix+1:end,:);
 estdates =estdates(ix+1:end);
 datestr(estdates(1,:))
 
 jmppoint = mod(3-month(estdates(1)),3); %to make sure that the first obs is end of quarter. 
 

 qestdata3 = estdata((jmppoint+1):3:end,:);
 qestdates =  estdates((jmppoint+1):3:end,:);

if mod( month( qestdates(1,:)),3) ~= 0;
    disp('Got the start date wrong')
end


stpoint = find( qestdates == datenum(2001,3,1));
enpoint = length(qestdates);
slevelvvalues= NaN(enpoint,5);
 
for zx = stpoint:enpoint;
    newsv = [0.0032  0.7078    0.1049    0.1328 0.1];
    rvll =   @(params)llf_restat(qestdata3(1:zx,:)' ,params,[],[],4);
    outvalues = fminunc(rvll,newsv , optz);
    slevelvvalues(zx,:) =    outvalues;
    
  
 
    
end




figure
subplot(3,2,1)
plot(qestdates,slevelvvalues(:,1))
title('Risk Premium Persistence')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,2)
plot(qestdates,slevelvvalues(:,2))
 title('Temp Shock Persistence')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,3)
plot(qestdates,abs(slevelvvalues(:,3)))
title('Perm Shock Std. Dev')
 datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,4)
plot(qestdates,abs(slevelvvalues(:,4)))
title('Temp Shock Std. Dev')
 datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])


subplot(3,2,5)
plot(qestdates,abs(slevelvvalues(:,5)))
 title('Risk Prem Shock Std. Dev')
datetick('x')
vaxis = axis
axis([qestdates(stpoint) qestdates(enpoint)+30 vaxis(3:4)])




orient tall
print -dpdf simplestuffinLEVELS


